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PHAZR:  A  Hienonenological  Code  for  Holeboring  in  Air 
I.  Introduction 

This  report  describes  a  new  code  for  studying  holeboring  by  a  charged 
particle  beam,  laser,  or  electric  discharge  in  a  gas.  The  coordinates  which 
parcdieterize  the  channel  are  radial  displacement  (r)  frcm  the  channel  axis 
and  distance  (z)  along  the  channel  awis  from  the  energy  source,  hence  the 
letters  "ZR"  in  the  name.  The  code  is  primarily  "phencmenological";  that 
is,  we  use  closed  solutions  of  simple  models  in  order  to  represent  mary  of 
the  effects  which  are  important  in  holeboring.  In  the  case  of  a  particle 
beam  propagating  through  a  gas,  we  use  a  model  provided  by  Lee  and  Buchanan 
[17]  to  account  for  changes  in  the  area  of  energy  deposition  (beam  radius) 
with  distance  from  the  accelerator.  We  have  modified  the  treatment  of 
scattering  of  the  beam  by  the  ambient  mediun  to  include  the  results  of 
Hughes  and  Godfrey  [18] ,  as  shown  in  i^pendix  B,  lb  represent  hydrodynamic 
expansion  of  the  gas  heated  by  a  given  energy  pulse,  we  use  the  adiabatic 
equation  of  state,  which  Latipe  et  al.  have  shown  to  be  a  good  approximation 
[1-3] .  The  equations  of  Picone  and  Boris  [4]  form  the  basis  of  a  subgrid 
turbulence  model  required  to  compute  the  enhanced  channel  cooling  caused  by 
asymmetry-driven  turbulence.  The  numerical  simplicity  which  we  gain  from 
the  use  of  such  solutions  enables  us  to  estimate  the  structure  of  a  channel 
while  using  a  far  less  computer  time  than  a  more  "detailed"  code  which 
solves  differential  equations  for  the  fields,  beam  envelope,  and  the 
chemistry  and  dynamics  of  the  ambient  gas.  This  feature  permits  the 
confutation  of  channel  properties  over  long  propagation  distances  (e.g, , 
the  distance  over  which  significant  beam  expansion  occurs)  and  thus  makes 
PHAZR  a  useful  code  for  those  studying  and  designing  future  systems. 

Manuscript  approved  June  13, 1986. 


As  indicated  by  the  report  title,  we  have  been  interested  primarily  in 
air,  although  the  model  will  apply  to  ary  chanically  nonreactive  gas  and  can 
be  modified  to  include  the  affects  of  chemical  reactions.  Conversely,  the 
various  modules  conprising  the  code  can  be  incorporated  easily  into  other 
codes.  To  account  for  tiie  effects  of  air  chemistry  in  PHAZR,  we  currently 
use  a  real  air  equation  of  state  routine  based  on  the  equilibrium  data  of 
Gilmore  [5-6] ,  vrfiich  cover  a  temperature  range  of  at  least  300  -  24000  K. 
Using  a  fast  table  lookup  routine  developed  by  Young  [7] ,  we  are  able  to 
obtain  values  of  y,  the  ratio  of  principal  specific  heats,  over  the  entire 
grid  in  roughly  the  same  amount  of  time  required  to  compute  a  vector  square 
root  function  on  the  Texas  Instruments  Advanced  Scientific  Oonputer  at  NRL. 

From  this  discussion,  we  see  that  the  name  PHAZR  forms  an  acronym  for 
"Phencroenological  Holeboring  in  Air  using  (Z,  R)  coordinates".  In  the 
remaining  sections,  we  discuss  the  various  modules  in  the  code  in  greater 
detail.  The  primary  emphasis  of  this  report  will  be  charged  particle  beams, 
and  as  an  exanple,  we  present  typical  results  for  an  ETA-1 ike  beam 
propagating  in  air.  These  calculations  will  clearly  demonstrate  how  PHAZR 
may  be  used  to  investigate  accelerator  parameter  space  and  to  isolate  the 
inportant  physical  parameters  which  determine  the  holeboring  properties  of  a 
given  system, 

Vfe  note  that  the  subgrid  tuirbulence  model  is  of  particular  interest 
because  of  two  features: 

(1)  This  model  represents  our  first  attempt  at  a  systematic  application 
of  the  analysis  of  asynmetry-generated  turbulence  to  a  numerical  model  of 
channel  physics. 


(2)  We  provide  a  link  between  the  theory  of  Boris  and  Picone,  which 
deals  with  the  generation  of  large  scale  turbulent  structure  (vorticity)  and 
the  e^qjerimental  data  of  Greig  et  al.  [11] ,  in  which  the  cascade  to  snaall 
scales  has  taken  place. 

The  problem  of  multipulse  energy  deposition  leads  naturally  to  a 
representation  in  which  the  effective  turbulent  diffusivity  varies  with 
position  in  the  interior  of  the  channel,  since  the  pulses  most  likely  will 
not  be  collinear.  Certainly  more  turbulent  transport  will  occur  near  the 
paths  of  the  pulses  than  in  portions  of  the  channel  which  do  not  contain  the 
trajectory  of  at  least  one  pulse.  In  addition  fluid  dynamic  and  turbulent 
transport  will  affect  the  distribution  of  turbulence  amd  the  spatial 
variation  of  the  turbulent  diffusivity  within  the  channel.  The  model  in 
this  report  has  the  above  properties  and  should  provide  a  useful  starting 
point  for  future  treatments  of  asyninetry  -  generated  turbulence. 
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II.  General  Structure  of  the  Code 

Because  we  assune  that  a  given  beam,  discharge,  or  laser  channel  is 
locally  cylindrically  symmetric,  we  use  only  the  radial  (r)  and  axial  (z) 
coordinates.  For  this  discussion  we  assume  that  we  have  a  model  with  closed 
solutions  or  tabulated  data  to  describe  energy  deposition  as  a  function  of 
axial  position  (z)  and  time,  i^pendix  B  gives  an  example  of  such  a  model 
[17]  for  a  particle  beam.  The  most  inportant  beam  parameters  for  energy 
deposition  at  a  given  z  are  pulse  radius  and  length  and  average  particle 
energy,  all  of  which  we  may  estimate  fron  a  model  or  data,  which  depend 
mainly  on  local  ambient  conditions  at  the  time  a  pulse  propagates  past  a 
given  point  cind  not  on  subsequent  channel  evolution.  In  order  to  compute 
the  local  conditions  for  each  beam  pulse,  we  need  to  perform  calculations  of 
channel  expansion  only  on  transverse  slices  of  the  channel  which  are 
displaced  in  z  by  distances  that  are  small  ccropared  to  the  shortest  scale  of 
variation  of  the  beam  parameters  given  above.  For  example,  for  a  particle 
beam  in  a  gas,  the  Nordsieck  length  normally  defines  the  shortest  scale  on 
which  the  beam  pulse  expands,  and  we  would  choose  our  slices  at  intervals  of 
seme  fraction  of  the  Nordsieck  length.  The  code,  therefore,  is  something 
less  than  two  dimensional,  since  the  numerical  grid  consists  of  a  nuttber  of 
weakly  connected  radial  grids,  each  at  a  different  value  of  z. 

For  a  given  slice  (transverse  plane)  of  the  channel  and  a  given  pulse, 
the  calculation  proceeds  as  shown  in  Fig.  1.  In  this  paragraph,  we  will 
briefly  discuss  each  box.  The  numbers  correspond  to  those  in  Fig.  1. 
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(1)  We  new  consider  a  string  of  N  pulses  and  a  set  of  J  "channel 
slices"  or  planes  transverse  to  the  path  of  propagation,  along  which 
energy  is  deposited  in  the  ambient  median.  We  compute  the  channel 
properties  consecutively,  beginning  with  pulse  1  and  slice  1,  after  the 
ambient  values  have  been  entered  during  initialization.  For  a  given 
pulse  "n",  we  proceed  to  each  channel  slice  "j",  canputing  the  evolution 
of  the  channel  up  to  the  time  at  which  another  pulse  is  to  be  produced 
fcy  the  accelerator. 

(2)  When  we  begin  a  new  pulse  n,  the  simulated  time  elapsed  since  the 
passage  of  the  first  pulse  through  slice  j  is  (n-1)  times  the  pulse 
repetition  interval  (PRI).  Note  that  the  time  step  nunber  "i"  may  be 
different  for  each  slice,  since  the  time  step  size  is  determined  by 
local  conditions, 

(3)  Fetch  the  values  of  the  necessary  variables  (density,  grid  cell 
positions,  tenperature,  ratio  of  specific  heats,  and  vojrticity)  from  memory 
for  slice  j. 

(4)  From  the  vorticity  ,  we  compute  the  turbulent  drffusivity 

and  then  a  turbulent  thermal  conduction  coefficient  X_. .  We  then  edd 

Ti  Ti 

to  the  classical  moleculcir  thermal  conductivity  based  on  kinetic  theory 
[8,9]  to  obtain  a  total  effective  thermal  conduction  coefficient  X^,  From 
the  temperature  and  the  ratio  of  specific  heats,  we  compute  the  energy 
density  at  each  radial  grid  point, 

(5)  Modify  the  vorticity  distribution  based  on  the  turbulent 
diffusivity 
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(6)  Calculate  a  time  step  size  based  on  the  timescale  of  energy  transport 
by  the  effective  thermal  conduction  corputed  in  (4). 

(7)  Vie  test  to  see  if  we  are  just  beginning  the  calculation  on  slice  j 
for  pulse  n. 

(8)  If  so,  we  ccnpute  the  radial  distribution  of  energy  deposition  by 
pulse  n  arri  add  this  incranent  to  the  existing  energy  distribution. 

In  the  case  of  a  particle  or  laser  beam,  this  step  includes  the  beam 
dynamics,  which  determines  the  radius  of  energy  deposition  as  a  function 
of  propagation  distance  (slice  nunber) .  This  factor  provides  a  physical 
link  between  successive  slices  along  the  propagation  path. 

(9)  If  the  test  in  (7)  is  positive,  we  also  ccnpute  the  increment  in 
vorticity  caused  by  pulse  n,  but  we  do  not  yet  update  the  vorticity 
distribution. 

(10)  We  modify  the  energy  distribution  based  on  the  effective  thermal 
conduction  which  we  calculated  in  step  (4). 

(11)  The  thermal  conduction  process  and  the  periodic  energy  deposition 
by  bean  pulses,  laser  pulses,  or  discharges,  cause  a  departure  frcn 
pressure  equilibriun  in  this  (time-step-split)  algoritlm.  Wfe  return  to 
pressure  equilibriun  in  one  time  step  by  using  the  adiabatic  equation  of 

state,  P  =  A  p  ,  where  P  is  the  pressure,  p  is  the  mass  density,  t  is  the 

ratio  of  principal  specific  heats  and  A  is  a  constant.  This  equilibration 
results  in  the  transporting  of  mass,  energy,  and  vorticity  in  a  Lagrangian 

manner.  This  transport  will  be  largest  when  a  new  pulse  propagates  through 

the  slice  (i.e. ,  t  =  (n-1)  PRI). 

(12) -(13)  If  this  is  the  first  time  step  for  pulse  n  (t  =  (n-1)  PRI), 


add  the  incrannent  in  vorticity  conputed  in  (9)  to  the  total  vorticity. 

Vte  do  this  after  the  fluid  transport  step  (11)  above,  since  the  vorticity 
increraent  which  we  ccnpute  is  the  value  present  after  the  hot  gas  has 
reached  pressure  equilibrium. 

(14)  L^ate  the  time. 

(15)  Perform  diagnostics. 

(16)  -(18)  When  a  time  inten/al  equal  to  the  PRI  has  passed,  store  all 
of  the  primary  variables  (see  (3))  and  increase  the  slice  index  by  one. 
(19)-(2C)  When  all  J  slices  have  existed  for  a  relative  time  interval 
n  X  PRI,  we  begin  the  sequ'J'  ce  for  a  new  pulse,  starting  at  slice  1, 
which  is  nearest  the  accelerator. 

(21)  When  all  pulses  have  propagated,  stop. 

In  the  remaining  sections,  we  present  the  equations  underlying  the  sequence 
of  steps  listed  above.  The  groupings  around  which  the  paper  is  organized 
are  as  follows:  energy  deposition,  thermal  conduction,  fluid  transport,  and 
subgrid  turbulence.  The  section  on  thermal  conduction  includes  a  discussion 
of  the  time  step  calculation,  and  that  on  subgrid  turbulence  includes 
various  aspects  of  vorticity  generation  and  evolution. 


III.  Deposition  of  Energy 

1.  Energy  Density  vs.  Radial  Distance 


Past  holeboring  calculations  have  employed  a  circuit  model  (Appendix  A) 
for  energy  deposition.  An  important  problem  with  circuit  nxDdels  relates  to 
the  radial  profile  of  energy  deposition.  Although  Bennett  profiles  with 
different  characteristic  radii  for  direct  and  ohmic  contribution  will  often 
provide  a  useful  representation,  detailed  models  like  VIPER  [15]  are  more 
realistic.  Ehirthermore,  we  must  account  for  the  decrease  in  particle  energy 
and  the  increase  in  beam  radius  with  increasing  propagation  distance,  as 
'A/ell  as  the  variation  of  energy  deposition  and  net  current  with  the  local 
density  at  a  given  position  along  the  beam  path.  Constructing  and  using  a 
table  of  ohmic  and  direct  energies  versus  these  parameters  would  entail  too 
much  expense  and  operational  difficulty  relative  to  the  accuracy  which  we 
would  expect  frcm  the  overall  model. 

We  therefore  combine  the  best  features  of  ’’detailed"  models  and 
circuit  models  by  using  the  following  schane: 

(1)  For  a  short  distance  from  the  accelerator  (~1  m) ,  use  a  "first 
principles"  code  like  VIPER  to  compute  tables  of  direct  and  ohmic  energy 
density  vs.  radial  displacement  from  the  beam  axis  vs.  local  mass  density. 

(2)  At  farther  distances  z  >  1  m  along  the  beam  path,  use  the  local 
average  mass  density  to  ccrpute  a  corresponding  net  current  from  the  Viper 
data  and  then  compute  the  average  beam  radius  a(z),  pulse  length  L(z),  and 
particle  energy  from  a  beam  dynamics  model  or  data. 


(3)  Using  the  local  density  at  z  >  1  m  and  the  VIPER  data  correspond j,r^ 
to  a  1  m  propagation  distance,  calculate  (interpolated)  radial  distributions 


of  direct  and  ohmic  energy  densities,  Oc^y/av)  and  OS^/av),  \(rfiere  d^^ 
denotes  the  amount  of  energy  deposited  through  direct  collisions  within  the 
volume  element  dv  and  d?  carries  a  similar  definition  for  ohmic  heating. 

(4)  Now  assume  that  the  energy  deposition  densities  at  axial  position 
z  >  1  m  have  the  forms  fj^(r/a(z)),  and  f^(r/a(z)).  Wfe  may  then  preserve  the 
shapes  of  the  functions  confuted  in  step  (3)  above  while  accounting  for  beam 
expansion  and  conserving  energy  throt^h  the  transformation 
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(III.l) 


In  eg.  (III.l),  the  "1"  signifies  the  first  slice  (at  approximately  1  m) , 
for  which  we  actually  have  data,  and  r'  is  the  radial  coordinate  for  the 
channel  slice  at  z.  This  gives  us  the  total  direct  energy  deposition  per 
unit  length: 

3^  (z)  =  27r  /  ^  (r',  z)  r'dr'  =  2Tr  /  (r,  1)  r  dr  (III. 2) 


We  use  a  similar  transformation  for  the  density  of  energy  deposited  at 
z  >  1  m  through  ohmic  heating.  Notice  that  we  have  not  yet  accounted  for 
the  decrease  in  pulse  length  with  increasing  z. 

(5)  Scale  the  direct  energy  density  to  account  for  the  decrease  in 
particle  energy  and  pulse  length  with  increasing  propagation  distance 
according  to  eq.  (B13),  developed  in  ^pendix  B  from  the  circuit  model 


(Appendix  A).  This  gives  us 


as 


as 


dv 

where 


ST  (|rff '  dV  1'  t-^l  [^1  1)  (III.3) 

a-^(z) 


R(z,  1)  =-| 


TiTTrr 


(III.4) 


is  the  ratio  of  the  respective  changes  in  particle  energy  per  unit  path 
length  due  to  direct  deposition  (subscript  D)  in  ambient  air  (p  =  p^)  for 
particles  of  energy  £(z)  and  e(l).  in  eq.  (III. 3)  L(z)  is  the  pulse  length 
at  propagation  distance  z  fron  the  accelerator. 

(6)  Again  using  the  circuit  model,  scale  the  ohmic  energy  density  dis¬ 
tribution  for  the  effective  inductance  according  to  eqs.  (B16)  and  (B17). 
This  gives  us 


dv 

where 


^  ^  r,  z)  ^  (r,  1)  [4511]  S(z,  1)  ( 

a^  (z) 


III. 5) 


In  (b/a(z)) 

Tn  (Vad)}  '  ^ 


(III. 6) 


S(z,  1)  = 


,  a(z)  >  b 


Here  b  is  the  radial  distance  frcm  the  axis  of  the  beam  at  which  the  conduc¬ 
tivity  falls  to  a  negligible  value  and  is  assuned  to  be  constant.  This  will 
be  a  useful  approximation  if  the  time  dependence  of  the  induction  reduces 
the  dmic  energy  deposition  by  a  constant  factor  o,  as  in  eq.  (AlO). 


2.  Energy  Deposition  Off  Axis 

lb  discuss  energy  deposition  by  a  pulse  having  a  trajectory  vdiich  lies 
off  the  channel  axis,  we  first  must  recall  that  we  have  tables  or  equations 
for  energy  density  (denoted,  in  this  section  only,  by  5)  deposited  versus 
radial  displacanent  (r)  fran  the  center  of  the  pulse  (Rj^)  versus  the  local 
mass  density  (p)  in  the  region  where  energy  is  being  deposited.  The  usual 
procedure  is  then  to  use  the  actual  local  density  (and  if  tables  are  being 
used)  to  interpolate  linearly  in  density  to  obtain  c(p)  versus  (r  -  Rq). 
When  ^  0,  a  one-dimensional  cylindrical  model  can  only  treat  energy 

deposition  as  occurring  in  an  annulus  centered  at  Rq.  Thus  we  initially 
treat  the  table  as  data  on  c(p)  versus  (r  -  Rq),  that  is,  with  r  and  Rq 
scalar,  to  interpolate  onto  our  computational  grid. 

Note,  however,  that  an  cinnulus  of  characteristic  width  "a"  at  sane 
radius  Rg  *  0  will  have  a  greater  area  than  a  circle  of  radius  a,  on  which 
the  tables  are  based.  We  must,  therefore,  renormalize  the  energy  deposited 
on  the  grid  by  the  factor 


[q  5(r,  0,  a)  r  dr 


(III. 7) 


^  [q  C(r,  Rq  ,a)  r  dr 

where  r  is  the  radial  coordinate,  Rq  is  the  radial  coordinate  of  the  center 
of  the  latest  pulse,  and  a  is  the  Bennett  radius  of  the  pulse.  For  exaitple. 


an  annular  Bennett-like  profile 


f(r,  Rq,  a)  = 


(III. 8) 


gives  us 


p  p 

R^  (a,  Rq)  =  { 1  +  “  f  +  tan-l  (-|)  ] }“! , 


(III. 9) 


vv 


which  is  always  <  1.  If  the  grid  is  sufficiently  extensive,  the  denoninator 
in  eq.  (I II. 7)  is  approximately  equal  to  the  integral  of  the  unrenormalized 
energy  density  over  the  grid.  If  not,  we  assume  that  the  energy  deposition 
has  the  form  of  eq.  (III. 8)  and  use  the  value  given  by  eq.  (III. 9). 
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IV.  Therroal  Conduction 

The  only  time  dependent  physical  process  currently  in  the  code  is 
thermal  conduction,  since  we  are  primarily  interested  in  the  evolution  of 
the  channel  temperature  after  expansion  to  pressure  equilibrium.  The  design 
of  PHAZR  assumes  that  the  time  scale  of  the  latter  process  is  short  compared 
to  the  pulse  repetition  interval  and  the  time  scale  of  thermal  conduction, 
we  use  the  adiabatic  equation  of  state  to  set  the  channel  density  to  the 
hydrodynamically  equilibrated  value  inmediately  after  one  time  step.  Recent 
calculations  using  the  HINT  code  [141  have  shown  that,  with  a  nonequilibrium 
chemistry  and  a  time  dependent  fluid  transport  algorithm,  the  minimum  pulse 
repetition  interval  for  which  pressure  equilibrium  is  reached  at  the  end  of 
each  pulse  is  approximately  50  us.  This  value  is  also  the  effective  minimum 
for  which  turbulence  would  have  any  effect  on  an  ETA  or  ATA  beam  of  ten 
pulses.  Thus  PHAZR  appears  to  be  useful  for  estimating  the  effects  of 
turbulence  on  channel  temperature  in  most  experiments. 

Cur  thermal  conductivity  is  the  sun  of  two  terms:  a  coefficient  of 
molecular  thermal  conduction,  computed  fron  kinetic  theory,  and  an  effective 
coefficient  of  turbulent  thermal  conduction  computed  as  described  in 


section  VI.  The  coefficient  of  moleculaa:  thermal  conduction  is  [8] 


where 


m  e 


3  1/2 

xO  =  (S  [erg/(cm  s  K)] 

SI  (T  )  ” 


(IV.l) 


(IV.  2) 


and  the  Eucken  coefficient  is  [9] 


C  =  0.115  +  0.354 
e  ^Y— 1 ' 


(IV.  3) 
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In  eq.  (IV.  2),  a  is  the  collision  diameter  associated  with  air;  n  is  a 

normalized  collision  integral;  M  is  the  molecular  weight  of  air;  T  is 

* 

temperature  in  K;  and  T  is  a  normalized  temperature.  In  eq.  (IV. 3),  y  is 

the  ratio  of  principal  specific  heats. 

Given  the  total  coefficient  of  thermal  conductivity  X  =  X  +  X_,  we 

m  T 

solve  the  diffusion  equation, 

If  =  7  •  X7T  (IV.4 

explicitly  using  cell-centered  differences,  where  5  is  the  internal  energy 
density.  The  choice  of  the  timestep  is  consistent  with  the  criterion  of 
Richtmyer  and  Morton  [10]  in  order  to  maintain  the  stability  of  the  scheme 
In  fact  we  use  the  somewhat  more  stringent  condition 

At  =  0.25  min  {N^k  ^r^/  -  1)}.  (IV.5 

In  eq.  (IV.5),  we  are  taking  the  miniraun  of  the  quantity  in  brackets  over 
the  grid,  where  is  the  particle  density  in  cell  i,  A^  is  the  area 
enclosed  by  cell  i,  Ar^^  is  the  width,  is  the  circunaference  of  a  circle 
passing  through  the  cell  center,  and  X^  is  the  thermal  conductivity.  The 
quantity  in  brackets  is  approximately  the  time  required  for  thermal 
conduction  to  smooth  out  a  temperature  gradient  across  a  given  cell. 


V.  Fluid  Transport  Algorithm 


As  stated  in  the  previous  section,  we  are  primarly  interested  in 
channel  evolution  after  pressure  equilibriun  has  been  reached  following  the 
passage  of  a  pulse  through  a  given  slice  of  the  channel.  In  particular, 
PHAZR  is  applicable  to  multipulse  beams  with  a  pulse  repetition  interval  of 
approximately  50  ys  or  greater.  Thus  we  need  only  use  adiabatic  expansion 
to  estimate  the  change  in  channel  properties  when  a  new  pulse  has  propagated 
through  the  channel.  Over  the  time  scale  of  channel  expansion  to  hydro- 
dynamic  pressure  equilibrium,  we  assume  that  turbulent  energy  transport  has 
very  little  effect.  That  is,  the  time  scales  of  energy  transport  by 
turbulence  and  molecular  thermal  conduction  are  much  longer.  Therefore,  we 
allow  hydrodynamic  equilibration  to  be  instantaneous.  The  vorticity 
distribution  changes  due  to  this  hydrodynamic  transport,  as  well  as  energy 
and  mass  density.  After  pressure  equilibration  occurs,  only  the  turbulent 
velocity  field  remains.  Since  we  use  an  effective  diffusivity  (thermal 
conductivity)  to  model  the  turbulence  field  and  because  the  internal  energy 
density  is  much  greater  them  the  kinetic  energy  density,  we  may  eliminate 
the  calculation  of  radial  velocities.  The  expansion  of  the  channel  is 
Lagrangian,  giving  us 


Af  *  A?  (P./P.)^^"^  (V.l) 

where  P^^  is  the  pressure  prior  to  transport,  P^  is  the  ambient  pressure,  A^^ 
is  cell  area,  superscript  zero  (0)  indicates  the  value  just  before  transport 
and  superscript  f  indicates  the  value  just  after  transport.  The  pressure 
is  the  sun  of  the  ambient  pressure  and  the  incremental  changes  in  pressure 
resulting  from  energy  deposition  by  a  pulse  and  from  thermal  conduction. 


Notice  that  at  the  end  of  each  time  step,  the  pressure  is  always  equal  to 
the  ambient  value.  Past  calculations  with  flux  corrected  transport  have 
shown  that  deviations  from  ambient  pressure  are  a  few  percent  after  the 
shock  from  the  most  recent  pulse  has  decoupled  from  the  channel.  Tto  conpute 
the  new  internal  energy,  we  derive  a  value  of  gamma  from  a  table  of  data 
measured  by  Gilmore  [7,8]  and  ijse  the  equation 


C  =  P/(y  -  1) 


(V.2) 


The  mass  density  is  given  by  the  (constant)  mass  in  a  cell  divided  by  the 
new  cell  volume.  Finally  we  must  conpute  a  new  value  for  the  vorticity  in 
each  cell,  according  to  the  equation 


f  0  f/  0 

zi  zi  ■  r 


(V.3) 


Where  is  the  z-ccn^wnent  of  the  vorticity  vector  u  and  i  labels  the 
ccmputational  grid  cells.  This  vorticity  expansion  term  does  not  include 
the  effects  of  diffusive  processes,  which  are  also  present  in  the  form  of 
the  effective  thennal  conductivity.  We  discuss  diffusive  processes  in  the 


next  section 


VI.  Subgrid  Turbulence  Model 
1.  Experimental  and  Theoretical  Background 

The  experiments  of  Greig  et  al.  [11]  have  provided  strong  evidence  that 
any  asynmetries  in  energy  distribution  by  an  electric  discharge,  laser 
pulse,  or  particle  beam  propagating  in  a  gas  will  lead  to  enhanceraant  of 
two  or  more  orders  of  magnitude  in  the  cooling  rate  of  the  resulting 
channel.  Schlieren  photographs  of  laser  and  electric  discharge  channels 
[11]  indicate  that  the  rapid  cooling  is  related  to  turbulence,  and  as  the 
channels  cool,  self-similar  expansion  occurs  according  to  the  simple 
empirical  equation 

R^(t)  =  R^(t)  +  4o  (t  -  t)  (VI. 1) 

where  R(t)  is  the  radius  of  the  channel,  t  is  the  time  measured  frati  the 
instant  at  which  energy  deposition  begins,  t  is  the  time  at  ^ich  expansion 
to  pressure  equilibriun  is  conplete,  and  o  (assumed  to  be  spatially 
nonvarying  in  the  channel)  is  an  effective  thermal  diffusivity.  Experiments 
show  the  diffusivity  to  be  approximately  constant  for  a  significant  period 
of  time.  This  behavior  would  result  fnxn  the  random  walk  of  vortex 
filaments  whose  strengths  do  not  decay  rapidly. 

The  fact  that  an  effective  diffusivity  provides  a  useful  model  of  the 
turbulent  transport  of  energy  indicates  that  a  subgrid  turbulence  model 
should  be  sufficiently  accurate  for  nunerical  holeboring  codes.  In  addi¬ 
tion,  we  note  that  the  details  of  the  turbulent  structure  (spatial  distribu¬ 
tion  of  different  scale  lengths)  must  not  be  of  major  importance,  since  such 
information  would  probably  not  lead  to  a  spatially  constant  diffusivity. 

The  size  of  the  region  in  which  the  turbulence  is  generated  should,  however. 


appear  in  the  calculation.  We  point  out  that  this  scale  may  differ  from  the 
channel  size  and  that  the  turbulent  diffusion  may,  therefore  not  initially 
operate  evenly  throughout  the  channel  cross  section. 

A  realistic  subgrid  model  should  also  provide  for  the  evolution  with 
time  of  the  turbulent  flew  and  the  effective  turbulent  diffusivity  in 
accordance  with  the  equation  for  the  evoluton  of  the  vorticity  field, 

2 

=  a  7  0)  .  (VI. 2) 

dt  z 

Since  our  calculation  uses  time  step  splitting  to  treat  fluid  transport  and 
thermal  conduction,  we  satisfy  eq.  (VI. 2)  in  two  steps: 

(1)  On  the  fluid  transport  step,  treat  the  vorticity  as  a  conserved 
Lagrangian  variable,  satisfying  the  continuity  equation  just  as  the  mass 
density  does  (see  eq.  (V.3)]. 

(2)  Use  the  effective  total  diffusivity  o  to  diffuse  the  vorticity. 

We  may  then  use  the  updated  vorticity  to  define  an  updated  effective 
turbulent  diffusivity  and  an  updated  total  thermal  conductivity.  We 
describe  the  relationships  between  these  quantities  below. 

Picone  and  Boris  [4]  have  developed  a  detailed  theory  of  vorticity 
generation  by  asynmetric  energy  deposition  in  a  gaseous  medium.  The 
important  asymmetry  types  which  the  theory  treats  include  the  following: 

(1)  Nor.collinear  pulses  or  misaligrsnent  of  an  existing  channel  and  a 
given  pulse, 

(2)  A  pulse  with  a  noncircular  cross  section, 

(3)  Three  dimensional  asynmetries  leading  to  energy  deposition  along  a 


curved  axis,  and 


(4)  Fluctuations  or  nonun if ormitites  in  the  energy  contours  within  a 
pulse.  An  excellent  example  of  case  (4)  is  the  existence  of  hot  spots 
within  a  laser  pulse  due  to  the  presence  of  several  inodes  [11]. 

The  generation  of  vorticity  follows  the  equation 

doj  - 

-^+a»7»v  =  (jj«7v+(Vpx  7P)/p‘^  (VI. 3) 

where  v  is  the  fluid  velocity  and  P  is  the  pressure.  Any  deviation  frcra 
cylindrical  symmetry  will  lead  to  the  misalignment  of  the  gradients  in 
pressure  and  density  as  the  hot  channel  gas  expands  to  pressure  equilibrium 
with  the  ambient  gas.  The  source  term  in  eq.  (VI. 3)  will  be  nonzero, 
leading  to  the  formation  of  at  least  one  vortex  filament  pair,  although 
cases  (2)  -  (4)  will  lead  to  more  cotiplex  vorticity  distributions.  The 
strength  (or  circulation)  <  has  the  form 

•Cig  =  [R(t)  -  B(0))  In  (p^/0(t))  f.g  (VI.4) 

where  i  labels  the  asynmetry  class,  6  labels  the  vector  component,  Ujj,(~Cg) 
is  a  characteristic  velocity  of  expansion,  p^^  is  the  local  mass  density 
prior  to  energy  deposition  by  a  pulse,  p(t)  the  density  at  the  position  of 
the  center  of  the  pulse  after  pressure  equilibration,  and  f^g  is  a  form 
factor  usually  <  1.  For  the  two-dimens icMial  asymmetries  [classes  (1),  (2), 
and  (4)1  ,  the  form  factor  f..  is  non-negligible  only  for  8  =  z,  for  which 
the  vorticity  vector  is  parallel  to  the  channel  axis.  Although  the  form 
factor  is  supposedly  calculable  for  the  cases  of  practical  interest, 
experimental  pulses  and  discharges  most  likely  contain  seme  combination 
which  would  be  difficult  to  discern  and  measure.  For  this  reason,  we  use 
eq.  (VI.  3)  with  f..  varying  over  a  reasonable  range  of  values  and  calibrate 
the  model  with  detailed  two-dimensional  calculations  employing  FAST2D  [4] . 


These  FAST2D  calculations  involve  the  solution  of  the  inviscid  equations  for 
conservation  of  niass,  mcrnentun,  and  energy,  and  include  no  effective 
diffusivity.  Thus,  in  tvo  dimensions  we  resolve  the  actual  asyninetry- 
induced  large  scale  turbulent  structure  (vortex  filament  pairs). 
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2.  Iheoretical  Effective  Diffusivity  in  the  Region  of  a  Vortex  Pair 

We  must  now  relate  the  vortex  strength  to  the  effective  diffusivity. 

The  theory  of  Picone  and  Boris  shows  that  asymmetric  energy  deposition 
produces  one  or  more  vortex  pairs.  The  vortices  in  each  pair  have  a 
separation  given  approximately  by  the  scale  length  of  the  local  asymmetry 
which  generated  the  pair.  To  obtain  an  effective  diffusivity  in  the  region 
containing  a  given  vortex  pair,  we  assune  that  a  rapid  cascade  to  small 
scale  structure  occurs.  This  is  reasonable  given  the  experimental  data  and 
the  proven  viability  of  an  effective  constant  diffusivity. 

Conventional  wisdom  often  presented  in  conjunction  with  dimensional 
analysis  states  that  most  of  the  energy  resides  in  the  large  eddy  scales 
while  the  dissipation  occurs  at  the  smallest  eddy  scales.  We  infer  that  the 
flux  of  fluid  through  the  region  originally  containing  a  vortex  pair  will  be 
constant  after  the  cascade  to  smaller  scales  occurs.  The  smaller  eddies 
will,  however,  mix  the  cooler  ambient  gas  being  pulled  into  the  channel  with 
the  hot  channel  gas  much  more  effectively  than  would  the  original  vortex 
pair  driving  the  flow. 

To  derive  a  relationship  between  a  and  k,  we  consider  fig.  2,  in  which 
a  pair  of  oppositely  directed  vortex  filaments  of  radius  and  strength 
±  |ic|  are  contained  in  a  cylindrical  region  of  radius  S.  We  denote  the 
separation  of  the  vortex  centers  by  25.  The  y-axis  is  the  line  passing 
through  the  vortex  centers,  and  the  x-axis  lies  midway  between  the  vortices. 
The  fluid  velocity  along  the  x-axis  is  [4] 


The  rate  (per  unit  vortex  filament  length)  at  which  ambient  gas  flows  into 
the  region  of  radius  S  and  between  the  vortices  is 


yo  /  dy'^  V  (0,0)  ~  o_  Ik  1  2(5-R)/itS, 
^  -(5-R^)  ""  ^ 


(VI. 6) 


where  M  is  mass  and  we  have  assumed  that  the  flew  is  nonnegligible  only 
outside  the  radii  of  the  vortices.  Numerical  simulations  [4]  have  indicated 
that  the  value  of  is  approximately  6/2.  The  mass  entrainment  rate 
implied  by  eg.  (VI. 1)  is 


3^  R^(t)  -  R^(t) 

-  =  '^“Pa- 


(VI. 7) 


Equating  the  mass  entrainment  rate  of  eq.  (VI. 7)  with  the  large  scale  flow 


rate  of  eq.  (VI. 6)  and  setting  =  6/2,  we  obtain 


(VI. 8) 


for  the  theoretical  effective  diffusivity  a.  The  turbulent  contribution  to 
the  thermal  conductivity  is  then  [12] 


^T  = 


(VI. 9) 


where  c^  is  the  specific  heat  when  the  pressure  is  constant. 

Unfortunately  eq.  (VI. 8)  is  not  yet  adequate  for  our  purposes,  since 
the  effective  diffusivity  applies  only  in  the  region  of  the  vortex  pair  or, 
equivalently,  in  the  region  of  the  nonuniformity  which  generated  the  vortex 
pair.  Therefore,  we  must  localize  the  turbulent  transport  coefficient 
within  a  distance  S^(the  scale  of  the  region  containing  the  vortex  pair)  of 
the  position  R.  of  vortex  pair  i.  For  example,  we  may  use  a  step  function 


^  0  ,  |r  -  R.  1  >  S. 


('71.10) 


where  is  given  by  eg.  (VI.  8).  In  the  case  of  PHAZR,  we  have 
axisynimetry,  so  that  a  pulse  which  occurs  off  axis  appears  as  an  annulus. 

Similarly  the  effective  diffusivity  resulting  fran  a  pulse  off  axis  will  be 
localized  in  an  annular  region.  Ifc  do  this,  we  use  a  Bennett-like 
function, 


a^(r)  = 


r  -  R.  ,  - 

fi  +  ! 


(VI. 11) 


where  r  is  the  radial  cylindrical  coordinate,  R^  is  the  distance  of  pulse  i 
fran  the  the  origin,  and  is  the  scale  of  the  density  depression  resulting 
from  pulse  i.  In  eg.  (VI. 11)  we  have  used  the  factor  [see  eg.  (111.9)1 


s  Rp  (S.,  R.) 


(VI. 12) 


to  account  for  the  increased  area  covered  by  the  annulus  relative  to  the 
original  circular  region  containing  the  vortex  pair.  Note  that  the  average 
of  o.(r)  is  egual  to  the  effective  diffusivity  o.  in  eg.  (VI. 8). 


3.  Model  for  the  Turbulent  Field 


We  new  assume  that  the  turbulent  velocity  field  may  be  approximated  by 
that  of  a  superposition  of  vortex  pairs  having  a  distribution  of  scale 
lengths  and  strengths.  Thus  we  will  deal  primarily  with  the  vorticity 
field,  the  effects  of  which  are  being  modeled  by  the  turbulent  diffusivity. 
We  will  assume  that  pulse  "i"  has  just  propagated  through  a  slice  of  the 
channel,  generating  a  vortex  pair  of  strength  lic^^  I  given  by  eg.  (VI. 4).  The 
scale  of  the  region  containing  the  vortex  pair  is  S^,  which  may  be  estimated 
from  the  equation  for  adiabatic  expansion 

S?  =  a?(P./P,)^^^  ,  (VI. 13) 

X  i.  X  a 

where  a^  is  the  radius  of  energy  deposition  by  pulse  i  and  is  the 
pressure  due  to  rapid,  local  heating  of  the  channel  gas  by  pulse  i.  Eq. 

(VI. 10)  then  gives  the  contribution  of  pulse  i  to  the  turbulent  diffusivity 
present  in  the  entire  beam  channel.  We  convert  this  to  an  increment  of  the 
total  vorticity  by  multiplying  eq.  (VI. 11)  by  ^/S^.  This  factor  comes  from 
first  using  eq.  (VI. 8)  to  convert  diffusivity  to  vortex  strength  and  then 
dividing  vortex  strength  by  the  area  containing  the  vortex  pair  to  obtain 
vorticity  (units;  s"M .  The  new  total  vorticity  then  becomes 

0)  {r,t.)  =  a)(r,t.  -  At)  +  4  o.(r)/S?  (VI. 14) 

X  X  XX 

where  a)(r,tj^  -  At)  is  the  vorticity  in  a  slice  just  prior  to  the  arrival  of 
pulse  i.  The  average  diffusivity  within  the  beam  channel  is  then  given  by 
the  area  integral  of  u  divided  by  4Tr,  to  be  consistent  with  eq.  (VI.  8)  and 
the  picture  of  beam  channel  turbulence  as  the  superposition  of  vortex 


Bather  than  use  the  average  channel  diffusivity,  we  preserve  the 
spatial  dependence  of  the  vorticity  field  by  defining  the  effective  channel 
diffusivity  as  the  product  of  u  (r,t)  and  a  factor  <S>2/4.  Here  <S>  is  an 
average  of  the  scale  sizes  of  the  vortex  pairs  generated  by  the  beam 
pulses.  For  pulses  of  the  same  Bennett  radius,  've  expect  ~  ~  S  for 

any  two  pulses  i  and  j.  Thus  we  have 


a(r,t^)  -  S  u(r,t.)/4 


(VI. 15) 


This  is  consistent  with  the  limiting  two-dimensional  case  shown  in  fig,  3. 
Here  the  pulses  are  all  misaligned  and  sufficiently  distant  that  the  local 
vorticity  is  determined  by  the  distance  fron  the  nearest  vortex  pair.  Wfe 


u  (r,t)  -  To),  (r,  S. ,  t) 
£  1  ~  1 


where  we  use  the  step  function  of  eg,  (VI.  10) 


(VI. 16) 


4a/S?  ,  I  r  -  r.  I  <  S. 

“i  ^i'  t)  =  {  I  r  -  r  1  >0 

0  ,  ~ 


(VI. 17) 


In  eg.  (VI, 17),  is  the  diffusivity  in  the  region  of  pair  i,  as  given  by 
the  integral  over  the  area  of  that  region.  The  total  effective  diffusivity 


is  then 


^  «(r)  dA  =  S. ,  t)  .S?  =  ^  .(r,  t)  (VI.IS) 


if  the  spatial  scales  of  the  vortex  pairs  are  all  egual.  Eguation 
(VI. 15)  is  therefore  consistent  with  this  limiting  case. 


VII.  Numerical  Tests  and  Simulations 

We  have  performed  tv«)o-dimensional  nunerical  simulations  with  the  code 
FAST2D  [4]  to  calibrate  the  results  of  PHAZR.  In  both  codes,  we  used  a 
table  of  data  derived  froti  VIPER  for  energy  density  deposited  versus  local 
mass  density.  The  data  correspond  to  a  distance  of  approximately  1  m  from 
the  accelerator  nozzle,  Wfe  assune  that  the  beam  parameters  are  as  follows: 
(1)  10  kA  beam  current,  (2)  5  ntn  pulse  radius,  (3)  10  ns  pulse  length,  and 
(4)  50  MeV  beam  energy,  we  used  pulse  repetition  intervals  of  1  ms  and 
100  lis,  respectively,  and  strings  of  ten  pulses.  Due  to  the  expense  of  the 
two-dimensional  simulations,  we  computed  channel  properties  only  in  the 
first  channel  slice,  which  was  located  approximately  1  m  from  the 
accelerator  nozzle. 

The  most  important  difference  between  the  PHA2R  and  FAST2D  calculations 
is  the  treatment  of  the  rotational  flows  generated  by  asynmetric  energy 
deposition.  In  the  case  of  FAST2D,  the  noncol linearity  of  the  pulses 
generates  most  of  the  vorticity  and  the  cascade  to  smaller  scales  is 
restricted  by  the  grid  size  and  the  flux-corrected  transport  algorithn. 

PHAZR  relies  on  the  subgrid  turbulence  model  described  in  the  previous 
section,  and  the  form  factor  f  in  eq.  (VI. 4)  can  range  frati  0.0  (no 
turbulence)  to  approximately  1.0  for  the  cases  encountered  so  far  [11,13]. 

We  must,  therefore,  perform  several  calculations  with  PHAZR,  each  with  a 
different  value  of  f,  in  order  to  compare  the  two  models.  (We  point  out 
that  the  notation  in  the  figures  is  "F"  instead  of  "f".) 

Figure  4  shows  our  results  for  ten  pulses  with  a  pulse  repetition 
interval  of  100  us  and  a  range  of  values  for  the  form  factor  in  PHAZR. 

Notice  that  we  ran  PHAZR  calculations  for  strings  of  collinear  pulses  and 
noncollinear  pulses,  respectively.  Each  pulse  in  the  "noncollinear"  strings 


was  displaced  frati  the  origin  by  0.0,  0.5,  or  1.0  times  the  pulse  radius. 

The  FAST2D  calculations  considered  only  noncollinear  pulses,  since  in  that 
code  no  turbulence  would  be  generated  in  the  collinear  case.  Vfe  chose  the 
pulse  locations  so  that  a  clockwise  progression  occurred  around  the  first 
pulse  which  was  centered  at  the  origin.  Clearly  there  are  an  infinity  of 
patterns  frcm  which  to  choose. 

For  the  100  us  case,  we  see  that  the  noncollinearity  of  the  pulses  was 
more  important  for  keeping  the  channel  cool  than  was  turbulence  (f  ^  0).  In 
addition,  turbulence  had  no  apparent  effect  until  the  sixth  pulse,  v^iich 
occurred  at  t  =  500  us.  This  is  the  time  scale  of  the  large  scale 
rotational  motion  [4]  and  shows  that  the  subgrid  turbulence  model  retains 
this  feature  of  the  motion.  Finally  we  see  remarkable  agreement  between 
F.4ST2D  and  PHAZP  when  the  latter  uses  a  form  factor  of  0.5. 

In  fig.  5,  we  show  the  results  of  a  similar  calculation  for  a  pulse 
repetition  interval  of  1  ms.  Here  turtwlence,  when  compared  with 
noncollinearity,  has  a  larger  effect  in  cooling  the  channel  than  for  the 
other  case,  since  the  time  between  pulses  (energy  depositions)  is  a  factor 
of  ten  greater.  Once  again  the  PHAZR  calculation  with  f  =  0.5  gives  close 
agreement  with  FAST2D. 

Figure  6  shows  an  example  of  the  application  of  PHAZR.  Wfe  use  the  same 
parameters  as  above  with  a  pulse  repetition  interval  of  100  us  for  a  string 
of  40  pulses.  In  this  case,  we  have  slices  at  1,  10,  20,  and  30  m  frcm  the 
accelerator  and  a  form  factor  of  0.5.  The  time  required  for  the  channel  to 
achieve  a  given  temperature  increases  monotonically  with  distance  frcm  the 
accelerator.  We  also  see  that,  because  of  turbulence,  the  channel  does  seem 
to  achieve  a  somewhat  uniform  temperature  toward  the  end  of  the  pulse 


train 


VIII.  Sutmary 

For  a  particle  beam  accelerator  which  emits  pulses  into  a  gas,  PHAZR  is 
useful  for  estimating  channel  properties  over  long  distances.  PHAZR  is  also 
applicable  to  studies  of  the  cooling  of  laser  or  discharge  channels. 

Because  this  code  does  not  involve  the  integration  of  fluid  or  field 
equations,  the  running  time  is  very  short,  making  PHAZR  valuable  for 
"systans-level"  studies.  Calibration  of  the  turbulence  model  with  accurate 
two-dimensional  hydrodynamics  simulations  has  shewn  excellent  agreement,  and 
the  initial  calculations  have  shown  that  pulse-to-pulse  noncoil inear ity  can 
result  in  a  measureaibly  cooler  channel,  even  if  turbulence  is  not 
significant.  We  mention  that  the  use  of  an  equilibriun  chemistry  model  and 
instantaneous  deposition  could  lead  to  a  channel  which  is  too  cool,  even  if 
there  were  no  shocks  to  carry  energy  out  of  the  channel.  Ihe  trends  which 
the  model  indicates  should,  however,  be  quite  useful  in  assessing  channel 
properties. 
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Fig.  1.  Flow  chart  for  PHAZR.  An  expression  "A  (bold  arrow)  B"  in^ilies 

that  we  may  derive  B  from  A  or  "A  gives  us  B".  An  expression  "A 

(thin  arrow)  A'"  means  that  the  quantity  A  is  modified  in  value  to  A'  under 


the  influence  of  quantity  B.  TSie  numbers  1-21  beside  the  boxes  correspond 


to  the  numbers  in  section  II  of  the  text 


15 


DIAGNOSTICS 


in^slies  that  we  nay  derive  B  from  A  or  "A  gives  us  B".  An  expression 
(thin  arrow)  A'"  means  that  the  quantity  A  is  modified  in  value  to  A' 


B" 

"A 

under 


the  influence  of  quantity  B.  The  numbers  1-21  beside  the  boxes  correspond 
to  the  numbers  in  section  II  of  the  text. 


Fig.  2.  We  define  the  region  of  influence  of  a  vortex  filament  pair  hy  a 
circle  of  radius  S.  Ihe  center  of  the  circle  forms  the  origin  of  our 
coordinate  system.  The  quantity  defines  the  cross  section  of  each  vortex 
and  the  separation  of  the  filaments  is  26.  The  vorticity  vectors  of  the  two 
filaments  are  oppositely  directed  and  the  quantity  <  denotes  the  magnitude 
of  the  strength  of  each  filament. 


Fig.  3.  We  consider  a  channel  produced  1:^  several  pulses  which  are 
displaced  from  each  other  hy  a  large  distance  con5)ared  to  the  radii  ,  Sj 
etc.  In  this  exan^ile,  the  vorticity  is  constant  and  nonnegligible  only 
inside  the  radius  of  the  ith  vortex  for  all  values  of  E.  Since  the 
pulses  are  all  roughly  identical,  the  regions  of  nonzero  vorticity  are  all 
of  approximately  the  seune  size  S. 
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Fig.  U.  This  figure  shows  the  results  of  PHAZR  and  FAST2D  calculations  of 
naxifflum  channel  temperature  Just  before  deposition  of  energy  ty  each 
pulse.  The  beam  parameters  are  given  above  the  diagram.  The  solid  lines 
give  PHAZR  plots  for  different  form  factors  and  collinearity  properties. 
The  dots  are  the  results  of  the  corresponding  FAST2D  calculation. 
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5.  This  figure  shows  a  set  of  calculations  which  are  similar  to  fig 
xcept  that  the  pulse  repetition  interval  is  1  ins  instead  of  100  ps. 


CHANNEL  TEMPERATURE  vs.  DISTANCE 


temperature  at  different  values  of  z,  the  position  along  the  beam  path 


APPENDIX  A 


A  SIMPLE  CTPCUIT  MODEL  FOR  ENERGY  DEPOSITIiDN  IN  THE 
AIMOSPHERE  BY  AN  ELECTRON  BEAM 


1.  Direct  Collisions 

The  average  "direct"  energy  loss  rate  per  electron,  per  cm  of 
propagation  distance  is 


=  £  r— 1 

'-dz^D  D  ''dz'^  I> 


(Al) 


v/here  z  is  the  distance  of  the  head  of  the  pulse  from  the  source,  p  is  the 

average  channel  density,  is  the  ambient  density,  and  the  change  in 

de 

particle  energy  due  to  direct  collisions  in  full  density  air  is 

2.5  keV/cm.  The  subscript  "Do"  on  the  right  hand  side  implies  direct 

de 

collisions  (D)  in  ambient,  full  density(»)  air.  Thus  sets  as  an 

effective  electric  field,  and  the  energy  in  joule/cm  deposited  at  z  by 
direct  deposition  is 


^2 

dz 


{A2) 


In  eg.  (A2),  is  the  total  energy  deposited  via  direct  collisions, 

Ij^  is  the  electron  beam  current  in  amperes  and  is  the  pulse  length  in 
seconds.  For  simplicity,  we  may  assune  that  the  radial  energy  distribution 
has  a  Bennett  profile: 


in  which  is  the  Bennett  radius  and  d^^/dV  denotes  the  density  of  energy 

deposited  through  direct  collisions.  The  Bennett  radius  a^,  for  a  given 

beam  emittance,  depends  on  the  net  current  whiah  may  be  determined  by  a 

calculation  using  the  VIPER  computer  code  115].  The  ratio  I  /L  varies 

n  o 

slowly  with  the  position  ^  (measured  frcm  the  front  of  the  pulse),  except  at 
the  expanded  head.  Constant  average  values  for  a^  and  therefore, 
provide  an  adequate  representation  for  the  circuit  model.  For  a  lOkA  beam 
current,  a  10ns  pulse  length,  and  ambient  density,  we  find  d^^dz  =<  0.25 
J/cm,  which  is  comparable  to  the  ohmic  deposition,  as  we  shall  see  below. 

2.  Ohmic  Deposition 

As  is  conmon  with  circuit  models  of  beam  energy  deposition  in  a  gas,  we 
assume  that  the  conductivity  a  is  high  enough  that  the  displacanent  current 
can  be  neglected  and  that  no  net  (beam  and  air  plasna)  charge  exists  out  to 
radius  b  »  beam  radius  aj^.  A  further  assumption  is  that  b  is  a  constant 
multiple  of  aj^(z=0).  Maxwell's  equations  (mks  units;  volts,  amperes, 
meters)  then  give  us 

7  X  E  =  -  (A4) 

~  3t 

or 

^  dr-  .  -  dr’  ,  (AS) 


since,  for  r  <  a  ,  E  (r)  is  approximately  independent  of  r  over  the  beam 

L/  Z 


cross  section  (variation  <  20%).  We  also  have 


1  3 


(Z  X  B)  =  r  B  =  uoJ„ 


(A6) 


virtiere  yo  =  4Tr  x  10"'  H/m  and  is  the  r^t  current  density  along  the 
axis  of  the  pulse  (z-axis).  Equation  (A6)  results  in 


Bg(r) 


Uo 

2^r 


{A7) 


where  I  (r)  is  the  net  (beam  +  plasma)  current  within  radius  r.  For  r>  a^ 
n  ~  D 

we  have 


Bg(r) 


(AS) 


in  which  is  the  total  net  current  and  is  independent  of  r.  Using 
eq.  (A8)  in  eq.  (A5),  we  find  that 


E 


—  * 

2it  ^n 


in 


(A9) 


Notice  that  the  factor  (u(]/2ir)in(b/ajj)  plays  the  role  of  an  inductance, 
which  is  assumed  to  be  constant  in  most  simple  circuit  models.  At  the 
radius  b,  the  conductivity,  by  definition,  becomes  too  small  for  the  plasma 
to  provide  space  charge  neutralization.  For  r  >  b,  therefore,  E^  drops  to 
zero,  an  assumption  which  was  used  implicitly  to  integrate  eq.  {A4).  A 
representative  value  for  b  is  usually  around  20  a^  (z=0).  Unfortunately 
calculations  using  VIPER  have  indicated  that  the  inductance  changes  with  5, 
leading  to  a  considerable  reduction  (factor  of  2  or  3)  in  ohmic  energy 
deposition.  A  crude  representation  of  this  would  be  to  multiply  eq.  (A9) 
by  a  Factor  a  =  1-f,  where  f  represents  the  fractional  reduction  in 


38 


electric  field  resulting  fran  an  inductance  which  varies  with  5.  We 
might  expect  o  to  vary  with  the  density  of  the  gas  (pre-existing  channel) 
into  which  the  latest  pulse  is  injected.  Thus  we  have 

|Ej  *  in(b/a^)  .  (AlC) 

The  energy  which  is  extracted  from  the  beam  (and  eventually  deposited  in  the 
air)  per  meter  of  path  traversed  is 

in  which  is  the  duration  of  the  pulse  in  seconds  and  the  integration 
variable  r'  is  the  time  since  the  pulse  head  traversed  a  particular  point  in 
space.  We  could  alternatively  integrate  over  5  =  ct  -  z,  which  gives 
position  in  the  pulse  relative  to  the  beam  head. 

In  keeping  with  our  model  of  the  beam  in  which  each  pulse  has  a 
constant  beam  current  of  duraction  we  neglect  the  following; 

(1)  erosion  of  pulse  length,  which  depends  on  the  propagation  range,  but  is 
typically  only  a  small  fraction  of  the  pulse  length; 

(2)  the  rise  of  at  the  head  (justified  after  the  rising  part,  typically 
-  0.3ns,  has  eroded  away); 

(3)  the  fall  of  at  the  pulse  tail,  which  results  in  relatively  little 
olmic  energy  deposition,  since  the  conductivity  a  is  high,  and  E  is  small 


Fran  eqs.  (AlO)  and  (All),  we  find  that  the  energy  deposited  through 
ohmic  heating  is 

^  ^  Ib  (A12) 

The  Bennett  profile,  eq.  (A3),  provides  a  less  faithful  representation  of 
ohnic  energy  deposition  than  of  direct  energy  deposition,  although  past 
holeboring  calculations  have  used  the  Bennett  profile  in  the  absence  of 
better  data.  Fran  Viper  calculations,  we  find  that  the  profile  for 
decreases  less  rapidly  with  radius  than  that  for  5  If  a  Bennett  profile  is 
used,  we  would  expect  that  the  Bennett  radius  a^  should  be  chosen  to  be 
greater  than  a^.  As  explained  in  Section  III  of  this  report,  we  circunvent 
this  problem  by  retaining  the  general  profile  of  the  VIPER  data  at  all 
slices  along  the  beam  path, 

we  would  expect  the  above  model  of  c^ic  deposition  to  suffer  somewhat 

from  its  simplicity.  In  the  beam  head,  the  conductivity  is  low  even  within 

the  beam  radius.  Electrostatic  neutralization  is  incomplete,  and  radial 

electrostatic  fields  and  displacement  currents  exist.  A  detailed  treatment 

[16] ,  however,  gives  the  same  result  for  ohmic  deposition  except  that  I^  is 

replaced  by  an  effective  current  I^,  which  includes  electrostatic  effects. 

At  t'=  t  (the  tail  of  the  beam),  I  =  I  ,  and  at  x"  =  0,  I  =  0,  so  that 
p  e  n  e 

the  integral  in  eqs.  (All)  and  (A12)  remains  the  same.  We  note  that  much  of 
the  energy  in  eq.  (A12)  is  extracted  from  the  nose  of  the  beam,  where  the 
radius  is  expanded,  and  electrostatic  effects  must  receive  consideration. 


Ihe  subsequent  deposition  of  this  energy  in  the  air  occurs  primarily  through 


E  .  J  in  the  main  body  of  the  beam  (or  even  after  the  beam  has  passed  by  a 
~  ~p 


given  location) ,  where  j  is  the  plasma  current  density.  In  the  body  of  a 

P 


pulse,  the  radius  is  approximately  constant,  E^  is  relatively  independent  of 


r,  and  J  has  a  profile  quite  similar  to  that  of  the  beam.  This  provides 
P 


sane  justification  for  approximating  the  profile  for  ohmic  deposition  by  a 


Bennett  profile.  Storage  of  energy  in  the  associated  magnetic  field  causes 


a  delay  between  extraction  of  energy  from  the  beam  and  deposition  in  the 


i 


fv 
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APPENDIX  B 


In  this  Appendix,  we  seek  to  transfonn  the  equations  of  Lee  and  Buchanan 
[17]  (belcw  'Aie  will  designate  this  reference  by  "LB")  into  fonns  whidn  are 
coirpatible  with  the  use  of  data  from  calculated  by  VIPER  rather  than  fron  a 
circuit  model.  Fran  this  stan<^int,  VIPER  provides  shapes  for  the  direct 

and  ohmic  energy  deposition  and  values  of  and  ,  the  respective 

dz  dz 

total  energies  deposited  per  unit  path  length. 


1 .  Particle  Energy 

Denoting  the  energy  per  particle  in  MeV  by  e ,  the  equation  of  LB  for  the 
change  in  particle  energy  with  propagation  distance  z  is 

i  =  -Tr-(3|)o’  '«> 


where  a  is  the  "radiation  length"  and  (de/dz)  is  the  change  in  particle 
energy  due  to  direct  deposition  (ionization  loss  rate).  Wb  may  also  define  an 
effective  ohmic  loss  per  particle  (de/dz )jj  by  dividing  the  ohraic  loss  by  the 
number  of  particles  in  a  pulse.  In  eq.  (Bl),  the  initial  particle  energy  is 
0.511  y  (MeV);  y  =  l/(l-6^)^^^;  6  =  v/c;  and  v  is  the  average  velocity  of  the 
beam  particles.  IB  treat  Ap  and  C^)q  ss  constants  in  integrating  eq.  (Bl). 
Assuning  that  these  quantities  do  vary  slowly  with  z,  we  suggest  that  simple 
averages  of  the  values  at  the  limits  of  integration  would  be  more  appropriate. 
If  we  integrate  eq.  (Bl)  between  positions  and  along  the  beam  path,  we 
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obtain 


'<^1+1’  *  ^i+l  -  ^i>/<XR>l,  l+ll 
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where 


and 


^^R^i,i+1  ^  ■*■  ^R^^i+1^^ 


"'^dz^D^  i,i+l  ^  ^  ^dz^D^^i+1^^ 


(B2) 


(B3) 


(B4) 


Vvhere  no  ambiguity  is  possible,  we  will  drop  the  subscripts  i  and  i+1  on  the 
average  (bracketed)  quantities  in  subsequent  equations,  and  we  will  denote 
quantities  evaluated  at  a  given  position  by  the  corresponding  subscript.  Ou: 
integrals  will  always  occur  between  positions  z^^amd 


de  d  *^^0 

We  will  new  express  terms  of  g—  and  g~,  the  total 


direct  and  ohmic  energies,  respectively,  deposited  per  unit  path  length.  If 


Ng  is  the  total  number  of  elections  in  a  beam  pulse,  we  have 


de  ^  1^  d? 

3z  °  N  3?  ' 
e 


(35) 


where 


N  =  n  wa^  L. 
e  e  D 


(B6) 


In  eq.  (B6),  L  is  the  pulse  length,  a  is  the  beam  radius,  and  n  is  the 

L/  6 

average  electron  density.  The  beam  current  is  then 


in  which  e  is  the  electron  charge  and  c  is  the  speed  of  light.  Substituting 


eq.  (B7)  into  eq.  (B6) ,  we  obtain 


N  — — 

e  ec 


de  ec  dc 


dz  LIj^  dz  ’ 


To  be  consistent  with  LB,  we  express  eq.  (B9)  in  terms  of  MeV/m,  assuming 


that  ^  is  given  in  erg/cm! 


^  ^  dc  _1  (1.6022  X  10~19C)(  2.9979  x  IQSm/s)  ( lOOan/m) 


dz  ”  3z 


LI.  (10 7erg/J) (1.6022  x  10"l9j/ev) (lO^eV/MeV) 


=  2.9979  X  10~3  dc  , 


(BIO) 


where  L  is  given  in  m  and  is  the  beam  current  measured  in  amperes  (A). 
Again  to  be  consistent  with  LB,  we  express  in  terms  of  the  value  in 


kiloamperes  (kA)  to  obtain 


2.9979xl0“«  dc  , _ 


(BID 


Equation  (Bll)  applies  to  both  olroic  and  direct  deposition  of  energy; 


however,  further  discussion  of  the  dependence  on  pulse  length  is  necessary  for 


;■ 


the  remaining  sections.  Both  direct  and  ohmic  deposition  per  electron  depend 
explicitly  on  L”^ .  However,  the  circuit  model  in  Appendix  A  indicates  that 

«  L  while  depends  on  L  mainly  through  the  logarithmic  variation  of 
effective  inductance  with  beam  radius.  Thus  is  approximately 

independent  of  L  while  varies  primarily  as  L“^ . 

At  positions  distant  from  the  accelerator  nozzle,  we  may  use  the  circuit 

model  (Appendix  A)  to  modify  VIPER  data,  which  give  and  versus  p  at 
positions  close  to  the  accelerator  nozzle.  If  the  first  position  (i=l) 


corresponds  to  the  VIPER  data,  we  have  from  eqs.  (Al)  and  (A2) 

^  rde-.  _ 
dz  *  id?' 


so  that,  for  position  i  >  1  and  a  given  local  channel  density  p. 


(812) 


We  define 


f— 1 

„  _  idz'l>  1 

il  ^ 

'■dz''l>  1 


(813) 


(814) 


for  which  approximate  values  may  be  obtained  fran  the  particle  energies, 
and  corresponding  to  positions  z^^  and  z^>  Zy  Fran  eqs.  (BID  - 
(814),  we  find  that 


.de.,  _  2.9979  x  10~6  „ 

iHz  'p^  I^  (kA)  i3z"^^  ^il' 


(815) 


which  is  indeed  independent  of  pulse  length,  as  we  would  expect  for  energy 
deposition  by  direct  collisions. 


For  ohmic  deposition  in  air  of  a  given  local  density  p,  we  have  eq.  (A12) 


II  -  in  (h/a^). 


Defining 


in(b/aj^^) 

®il  -  irTTBTi^' 


and  assuming  that  b  is  approximately  constant,  we  have 

2.9979  X  10“® 


/■de'i  _  X 

■dz^Qi"  L,  1.  (kAi 


1  *-b' 


(B16) 


(B17) 


(B18) 


In  the  section  B.3  we  shall  see  that  a^  varies  approximately  as  exp  [A(Lj^- 

Hl' 


Lj^)],  where  A  is  a  constant.  The  quantity  S.,,  therefore,  depends  only  weakly 


on  L,  as  we  would  expect. 


2.  Pulse  Length 


LB  assume  that  pulse  erosion  accounts  for  the  effects  of  ohmic  deposi¬ 
tion,  as  well  as  the  other  mechanisms  suntiarized  by  eq.  (Bl).  If  we  assume 
that  the  average  number  of  electrons  per  unit  length,  dN^/dL,  is  constant  in 
time,  we  have  in  the  formulation  of  IB 


d 

dz 


,  eL  r  cCJeN  ,  ~n” 


3z^ 


I_e 


{B19) 


in  which  A  is  the  inductance  obtained  from  a  circuit  model,  such  as  that  in 

i^pendix  A;  (kA)  is  the  net  current  (beam  current  +  plasma  current);  and 

the  Alfven  current  I^is  defined  as  1.7xl0i8Y(kA)  or,  equivalently,  e(MeV) 

times  34(kA/MeV)  from  LB.  As  in  eq.  (Bl),  IB  treat  X  ,  f^)_,  and  A  I  e/I. 

Iv  ^QZ'^D  n  Pi 

as  quantities  which  do  not  vary  with  L  but  which  do  vary  slowly  with  z. 
again  si^gest  replacing  them  with  averages  as  in  eqs.  (B2)-(B4).  An 
equivalent  equation  to  eq.  (B19)  is 


d 

dz 


sL 


D 


(B20) 


Coitparing  (B19)  and  (B20),  we  find  that 


As  noted  above,  the  left  hand  side  of  Bq.  (B21)  does  not  vary  strongly  with  L 

and  we  conclude  that  neither  does  the  right  hand  side.  Recalling  the 

discussion  at  the  end  of  section  B.l,  (eqs.  (B12  -  B18';)  '//e  see  that 

(^)  is  indeed  approxL^nately  proportional  to  ,  consistent  with  this 
^  a 

conclusion.  Tto  obtain  the  remaining  equations  which  correspond  to  (B21)  from 
those  given  by  LB,  we  may  therefore  replace  the  quantity  M^/3A  by 


^^^dz^  ^  ^  ^  H^dz^ai  4+1  4z^a(i+l)'^' 

X  f  X  •  L 


(B22) 


which  accounts  for  the  slow  z  dependence. 


Fran  eqs.  (Bl)  and  (B19),  tB  obtain 

dL  _  _  ^4  _  _  L  ^dg<. 
dz  34e  e  '■dz'h 


Using  eqs.  (Bl)  and  (B23),  LB  also  obtain  an  equation  for  dl/de. 


(B23) 

From 


that  expression,  we  obtain  the  pulse  length  in  eq.  (B24); 


3.  Pulse  Radius  (a) 


LB  give  the  IsJordsieck  equation  as 


ln(a/el  )  =  -Ar-  / 
d2  n 


(B25) 


where  =«  3700  (MeV  -  kA) .  v'Je  inay  make  the  substitution  defined  by  eq.  (B23) 


in  the  !^rdsieck  equation  to  obtain 
d  in  (aj^/el^) 


<X^>  <!„>= 


(B26) 


-  7400/(<Xj^>  <I^> 


Integrating  Eq.  (B25)  between  positions  and  we  have 


.  ^i'ni 


y^exp  {7400(L.-  L.^i)/(<V 


An  alternative  form  is  necessary  if  we  wish  to  use  the  expression  of 
Godfrey  and  Hughes  [18]  for  the  Nordsieck  length  L^.  Here  we  have 


L  -  7,14  vy  In  (1090  vy) 
N  [ln(2970vY)]2 


(B28) 


where  Budker's  parameter  is  defined  by 


1.7x10**  8 


(B29) 


in  which  I  is  the  Alfven  current  limit  and  I.  the  beam  current  in 
amperes  (A), 

From  eq.  (B  25),  we  identify 
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